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A three-dimensional fluid model which includes the dispersive effect of electron inertia 
is used to study the nonlinear macroscopic plasma dynamics of small scale discrete auroral arcs 
within the auroral acceleration zone and ionosphere. The motion of the Alfven wave source 
relative to the magnetospheric and ionospheric plasma forms an oblique Alfven wave which is 
reflected from the topside ionosphere by the negative density gradient. The superposition of the 
incident and reflected wave can be described by a steady state analytic solution of the model 
equations with the appropriate boundary conditions. This two-dimensional discrete auroral arc 
equilibrium provides a simple explanation of auroral acceleration associated with the parallel 
electric field. Three-dimensional fully nonlinear numerical simulations indicate that the 
equilibrium arc configuration evolves three-dimensionally through collisionless tearing and 
reconnection of the current layer. The interaction of the perturbed flow and the transverse 
magnetic field produces complex transverse structure that may be the origin of the folds and 
curls observed to be associated with small scale discrete arcs. Late time transverse electric field 
power spectra tend towards a universal Jr^/3 spectral form. 
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i. Introduction 


This paper presents a model of the three dimensional evolution of small scale discrete au- 
roral arcs originating from Alfven waves. The model which includes dispersion due to electron 
inertia, describes the region of the magnetosphere usually referred to as the acceleration zone, 
and in addition, includes the interaction of the Alfven waves with a collisional inhomogenous 
ionosphere. Haerendel, 1983, described the concept of oblique Alfven waves and formulated a 
model of auroral arcs in which he was able to explain some aspects of the transverse electric 
field structure. In the present paper, it is shown that nonzero electron mass combined with the 
oblique Alfven wave model leads to a natural explanation of the origin of a steady state parallel 
electric field as well as a mechanism for three-dimensional instability of the arc. 

The most important results emerging from the study of the proposed model may be 
summarized as follows. The oblique Alfven wave model of auroral arcs [j Haerendel , 1983] is 
extended to include finite electron inertia. Analytic and numerical solutions of this model are 
found that prove the existence of non-dissipative steady state wave-like solutions which have a 
complete two dimensional description of the electric and magnetic fields. This description 
includes a parallel electrostatic field which drives the field aligned current and accelerates the 
electrons. These east-west aligned, north-south drifting equilibria are consistent with the 
observed magnetospheric structure of discrete auroral arcs and offer a simple and self-consistent 
explanation of the origin of parallel electric fields. The two dimensional arc equilibria are three- 
dimensionally unstable to collisionless tearing and reconnection to produce interesting structure 
in the electric field and current density. 

It is very important to understand the physical scale of the phenomena that is being 
investigated in the present paper. The terminology that exists in the literature is confusing in this 
regard so an attempt will be made to place previous work in proper perspective to the present. 
Hallinan and Davis, 1970 refer to a class of discrete auroral phenomena as “small scale auroral 
arcs” which are also termed “breakup arcs” by Goertz, 1981. These breakup or small scale au- 
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roral arcs have a transverse scale of the order of a kilometer and sometimes have such associated 
features as folds and curls in the terminology of Hallinan and Davis. The spatial scale is related 
to the collisionless or electromagnetic skin depth due to electron density. Other auroral 
phenomena which are also termed discrete auroral arcs have a transverse scale (10 to 100 km) 
which in some models depends upon the conductivities of the generator and the ionosphere 
[Muira and Sato, 1980; Lyons, 1980; Chiu and Cornwall, 1980; Sonnerup, 1980; Lysak, 1985; 
Lotko et al, 1987]. The relationship between the various types of discrete arcs is not well under- 
stood, but it is often observed that the larger type of arcs form first and then break up by forming 
a system of arcs which are often described as draperies or curtains. A scenario for this breakup 
was suggested by Goertz, 1981. 

The present work is concerned with the former type of discrete arc and the terminology 
‘discrete auroral arc’ adopted in this paper refers to arcs having a transverse scale of the order of 
the collisionless skin depth. The acceleration zone is defined to be the region where c 2 /co 2 L^ 2 
°e B/n maximizes, where L ± characterizes a geomagnetic flux tube area. For small scale arcs, 
the acceleration zone is of particular importance since that is where wave dispersion and the par- 
allel electric field maximize. 

The two dimensional interaction of Alfven waves with the ionosphere has been 
considered by several authors [Goertz and Boswell, 1978; Lysak and Carlson, 1981; Lysak and 
Dum, 1983]. There is much useful information that can be obtained in considering arc evolution 
in two dimensions, but it is only in three dimensions that nonlinear behavior plays a significant 
role. In three dimensions, an east-west aligned arc can be unstable to several potential 
instabilities. Seyler, (1988) in a preliminary study showed that finite electron inertia leads to a 
rapidly growing collisionless tearing mode of the Alfven wave when it is bounded by a perfectly 
conducting static ionosphere and has a transverse scale comparable to the electromagnetic skin 
depth. Seyler found that the ionospheric boundary conditions are very important in determining 
the nature of the Alfven wave evolution. Specifically, for a boundary condition corresponding to 
zero parallel current density, an MHD shear flow instability resulted rather than the tearing 
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mode. Given the sensitivity of the evolution upon the ionospheric boundary conditions, a 
simulation study using a more realistic ionosphere model is needed. 

The primary objectives of this paper are to propose a realistic dynamical model to 
describe the evolution of Alfven waves in the ionosphere-magnetosphere environment and to 
present numerical simulation results which show the model self-consistently explains the 
structural morphology associated with discrete arcs and possibly energetic particle acceleration 
as well. 

The paper is organized as follows. Section 2 is concerned with the mathematical model, 
its derivation and physical properties. The development of the dynamical model proceeds from 
first principles, relying upon assumptions appropriate for the inner magnetosphere and the 
ionosphere. In Sec. 3 a stationary two dimensional solution of the equations is presented that 
represents a simple picture of a small scale non-dissipative discrete arc with the associated 
transverse electric field, parallel current density and parallel electric field which accelerates the 
electrons and supports the parallel current. The results of three dimensional simulations are 
presented in Sec. 4 show that stationary arcs can be unstable to a collisionless tearing mode 
which may be responsible for the observed transverse structuring in the form of folds and curls. 
A summary and a discussion of the most important results is given in Sec. 5. 

2. DYNAMICAL MODEL 


2.1 Derivation 

A dynamical model of the plasma which is valid on spatial scales between say one 
hundred meters and 100 kilometers transverse to the geomagnetic field and which extends from 
the E-region ionosphere to several Earth radii into the magnetosphere would be sufficiently 
complete to contain most features of the interaction of Alfven waves with the ionosphere. Such a 
model, however, is intractable insofar as one can compute numerical solutions to the complete 
model without a simplification of the physics and a limitation of the spatial scales involved. In 
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this section a model is derived which we believe contains some of the more important aspects of 
magnetosphere-ionosphere interaction and which is numerically tractable. The basic model 
consists of the following set of equations, which are a version of the warm quasineutral two- 
fluid equations. 


B,n- t- V -(mi) = 0 


mn(d, u + u • Vu) = — Vp + — J x B - mnv u 

V-o 


( 1 ) 

( 2 ) 


3,B = -V xE 
V x B = [ l 0 3 

E + u x B = p 0 X. 2 [d, J + V • (uj + Ju)l + — J x B — — Vp + ji 0 r| J 

ne 2 ne 


(3) 

(4) 

(5) 


where n is the plasma density, u is the fluid velocity, m is the ion mass, v is the ion-neutral 
collision frequency, T| is the magnetic diffusivity, p is the total pressure given by p=n(T e +Tj) , 
with T e = Tj = constant. The collisionless plasma skin depth is X = c/C0p, J is the current 
density, B is the total magnetic field (geomagnetic plus perturbed) and the electric field E has 
both an electrostatic and electromagnetic component and it is determined by the generalized 
Ohms law, Eq. (5). Equations (1) - (5) are an adequate description of the macroscopic plasma 
behavior on spatial scales much larger than the ion gyroradius. They do in fact contain 
considerably more physics than is essential to describe large scale auroral arc behavior in the 
magnetosphere and upper ionosphere (F-region). 

A simplified model may be derived from a systematic expansion in the small parameter 
e s 6)/Q- «1. For the region of interest - 1-2 R e the relevant nondimensional physical 
parameters are taken to scale in e as follows: 
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Lj_/L|| ~ e 8B/B 0 ~ e 8 uAa ~ e c/L_LCOp ~ 1 

8u/L x Qi~e pi/L ± -e pg/L^-e 2 vg/v A -e 

For now the collisional transport parameters v and rj will be taken to be zero which is the case 
for the magnetosphere. Later they will be restored when the ionosphere is considered. 

An important consequence of the ordering assumptions is that the magnetic field may be 
split up into a geomagnetic component and a perturbed transverse part due to field aligned 
currents. Accordingly write B as 

B = zB 0 + zxVy (6) 

where \|/ is the perturbed flux function that obeys 

VV = \L 0 Jz (7) 

Assume the parallel flow velocity is small u z - eluj. We then have Vj_- - 0(e 2 ), 

which may taken to be zero. Define the stream function <|> through u = zx V<j>. Henceforth let 
V = V ± and all z derivatives will be explicitly displayed as 3 Z . To first order in e the 
perpendicular part of the Ohms law becomes, 

E X = -£ 0 V<1> (8) 

For the z component the leading order terms are of order e 2 , 

2 2 

£ 2 = -zxV(}> Vy + X (3, + zx V<|>- V)V ( 9 ) 
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Equations (11) and (12) are two equations for the scalar fields <J> and i|/. They reduce to the 

2 

Strauss equations [Strauss, 1976], often called reduced magnetohyrodynamics, for k = 0. The 
collisionless skin depth introduces a length scale into the Strauss equations. Additional physical 
effects associated with this scale include collisionless tearing and reconnection and a parallel 
electric field. The discussion on the dynamical invariants of Eqs. (11) and (12) which follows in 
the next section should provide insight into the role of the collisionless skin depth. 

lit should be pointed out that the ordering assumptions invoked do not include the kinetic 
Alfven wave discussed by Hasegawa, 1976. This is because for the region of the magneto- 
sphere of interest <2R g , the hybrid ion gyroradius p- = (T e /mjQj 2 ) 1//2 , is considerably less 
than the skin depth and is therefore not important. The type of Alfven wave that is of interest 
here is more properly called the inertial Alfven wave. The terminology used in Seyler, 1988 is 
incorrect in this regard, since that paper dealt with the inertial Alfven wave and not the kinetic 
Alfven wave of Hasegawa. Although the terminology kinetic has been used to refer to both 
types of waves, it will not be used here. 
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2.2 Dynamical Invariants 

The magnetospheric plasma evolution equations given by Eqs. (11) and (12) possess two 
global constant constants of the motion for periodic or perfectly conducting boundary 
conditions. The total energy H is conserved and consists of three terms corresponding to 
transverse kinetic energy, the transverse magnetic energy and the parallel electron kinetic energy 
due to the current carrying electrons. The form of H is, 

// = J rf 3 Jc||V<t >| 2 + |Vv | 2 + x 2 |v 2 v | 2 } (13) 

There is another global invariant which is related to the cross energy of ideal magneto- 
hydrodynamics and is a generalization that includes parallel electron inertia. We denote it by T 
and it is. 


r = J d 3 x{ |v<j> ■ vtH + x. 2 v 2 \|iv 2 <i>} 


(14) 


There exists a local invariant which is a generalization of the Alfven frozen flux theorem. 


<D = J (B-X. 2 V 2 B) dS = J(\i/-X. 2 V 2 v)z dt 

S c 


(15) 


where S is an open surface and C is the contour bounding that surface. Then generalized flux O 
is locally conserved, that is it simply advected by the flow and thus obeys, 

(a, + zxV<t>-V)<D = 0 (16) 

This invariant shows that the magnetic flux through a surface moving with the fluid is not 
conserved for X * 0. Therefore reconnection and tearing of magnetic field lines can occur 
without resistivity. 
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2.3 Linear Properties 

Consider the dispersion relation for the propagation of a kinetic Alfven wave in a 
homogeneous equilibrium. We readily find, 


co 


2 


*z 2 Va 

1 + kl X 2 


(17) 


where k denotes the perpendicular wavenumber. This dispersion relation shows that small 
transverse scale waves propagate slower than larger scale waves. The consequences of this fact 
are important for the formation of multiple arc systems. For a spatially localized disturbance the 
transverse spatial Fourier spectrum of the disturbance is broad. If there is sufficient time for the 
Fourier components to separate as they propagate along the magnetic field, then the large scale 
components of the disturbance reach the ionosphere first, followed by short scale components. 
The effect appears as the breakup of a large scale auroral arc and the formation of a multiple arc 
system. This scenario for the breakup of arcs was first suggested by Goertz, 1981. 

Now consider the effect of dispersion on the breakup arc after it has formed. The 
breakup arc has formed an equilibrium that varies in x and is constant in y. This equilibrium has 
an E x B flow in the y direction as well as a perturbed magnetic field also along y, where the 
ratio of E x /B y is of the order of the Alfven speed v A . It is well known that the flow speed must 
exceed the transverse Alfven speed v^ = By/V4ttmn in order that the flow be unstable to the 
Kelvin-Helmoltz instability , [ Chandrasekar , 1961]. So that that the ideal (no dispersion) 
magnetohydrodynamic KH instability should not occur in this equilibrium with a large growth 
rate. However a collisionless tearing mode can be unstable, since dispersion allows for a 
slippage between the flow and the perturbed magnetic field as is shown by Eq. (15) and Eq. (16). 
For the case of no equilibrium flow, Seyler, 1988 gave an expression for the linear growth rate 
of the collisionless tearing mode, which is essentially a re-scaling of the classical result of Furth 
et al , (1963) , and is also given in the review article by White, 1986. The growth rate is. 
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( 18 ) 



where A' is the logarithmic derivative of the outer solution for the perturbed magnetic flux 
function, t a is the transverse Alfven transit time with respect to the magnetic field produced by 
field aligned currents and I is a dimensionless integral whose numerical value is about 1. The 
important thing to note is the strong dependence of the growth rate upon the skin depth. For 
typical magnetospheric auroral arc conditions kA. = 1, A'X « 1, and t a « 0.2 sec, gives y ~ 5 
sec -1 , which is of the order of the observed growth rate for the formation of auroral curls, 
[Hallinan and Davis, 1970]. 

2.4 Magnetospheric boundary conditions 

The magnetospheric boundary corresponds to a generator with an internal conductance 
G. This boundary condition is identical to that considered by Lysak, 1985, 

Xf + \i. 0 B 0 L G § = f{x,y,t) ( 19 ) 


where Eq is the generator conductivity and f(x,y,t) is a prescribed function. As Lysak has 
pointed out this boundary condition corresponds to a pure current generator for Eq — > 0 and a 
pure voltage generator for Eq and a nonreflecting generator for Eq' 1 = p o v A (l + k 2 X 2 ) ly2 . 

The form of the function f is taken to be a function of x with small traveling wave 
perturbations in the y-direction and is. 


f{x,y,t) = fo(x)\ 


1 + 


N 

^C„sin(ifc n y - cor)sin(jt cr x - cot) 


n=l 


( 20 ) 


Two forms of f 0 (x) will be considered. One of which is a Gaussian and the other for 
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which f G ' (x) is a Gaussian. These two forms of f Q (x) are intended to model a V shock and an S 
shock respectively in the terminology of Temerin et al, 1981. Specifically the function used is 
f D (x) = exp[-(x - 7t) 2 /A 2 ], where 0 < x < 2 ti. 

In the simulations to be discussed in Sec. 4, the magnetospheric boundary condition is 
chosen to be a non-reflecting generator. For substorm related events, the physical generator is 
far out in the magnetotail where the slingshot effect launches Alfven waves towards Earth. In 
the auroral acceleration zone, there is no generator per se, but only the flux of waves from the 
distant generator. So that by choosing a non-reflecting source of Alfven waves at about 2 R , 

V 

the flux of Alfven waves from the actual generator can be modeled without causing non-physical 
secondary reflections of waves incident on the magnetospheric boundary due to the physical 
reflections from the ionosphere. 


25 Simple Ionospheric Boundary Conditions 

Mallinckrodt and Carlson [1978] considered the reflection of Alfven waves from a con- 
ducting ionosphere. For a conducting ionosphere, the reflection coefficient is given by 


p _ Rj ~ Z A 

R i + Z A 


( 21 ) 


The limiting cases of this expression are: (1) short circuit T = -1: =><}> = 0 for Rj — » 0, 

(2) open circuit T = 1: => = 0 for Rj — » »o, and (3) nonreflecting T = 0: => <)> = 

v A (l+k2a2)l/2 v/Bo 

For the naive model of the ionosphere in which Rj = 1/Ep, a calculation shows for the 
auroral ionosphere it is generally true that Rj « Z A . In this case T ~ -1, and the ionosphere 
acts as a short circuit. This assumes, however, that Alfven wave reflections are due to the high 
conductivity of the ionosphere. As it turns out, a more important source of Alfven wave 
reflection is due to the negative density gradient and the associated increase in the Alfven speed 
encountered in the transition from the ionosphere to the magnetosphere, [Ellis and Southwood, 



1983]. This density gradient also gives a reflection coefficient of T ~ -1. The fundamental dif- 
ference between the good conductor boundary condition due to small Rj and the good reflector 
boundary condition due to a density gradient is that the good reflector boundary condition alone 
does not allow a cross field current to close the parallel current. A realistic model of the 
ionosphere should include a density gradient as well as ion-neutral collisions for current closure. 

3 . Steady State Auroral Arcs 

The model given by Eqs. (11) and (12) admits nontrivial two-dimensional steady state 
solutions if the the Alfven wave source is drifting normal to the plane of the current sheet and if 
certain magnetospheric and ionospheric boundary conditions are imposed. The stationary 
solution to be described is intended to represent a east-west aligned auroral arc during the 
expansion or recovery phase of an auroral substorm. During the expansion or recovery phase, 
the auroral oval is expanding or contracting in a north-south direction. This motion is 
represented by the drift velocity vj of the Alfven wave source. Let the geometry be as follows: 
x is north, y is east, and z is up. Transform to the frame of the drifting source, then x — > x - v d t. 
It is important to note that the plasma in the region of interest is not drifting with velocity v&, but 
only the source of Alfven waves. The drifting source of Alfven waves emits what Haerendel 
[1983] has appropriately called oblique Alfven waves. In steady state, d/dt = 0 and in two 
dimensions dy = 0, we have, 


vA4>+ " d;Y = o 

(22) 

v,a^-x J a, ! Y)+w i>=o 

(23) 


These equations may be combined to give 
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\j 8l(<j>-X. 2 d 2 <|>) = v A <J> 


(24) 


Before presenting the solution to Eq. (24), note that if X = 0, the form of the equation is a 
hyperbolic wave equation. Solutions to equations of this type are in general well posed if the 
boundary conditions are of open surface Cauchy type. For the problem at hand the physically 
relevant boundary conditions are closed surface Dirichlet, which in general is too restrictive for 
hyperbolic equations. The general solution of Eq. (24) without application of boundary 
conditions may be expressed in a Fourier series as 


<t>(*,z)= ^ 4)£expj/&r;t±5(l + X 2 & 2 ) 

k=- oo l L ' 


1/2 


(25) 


where ^ is the Fourier transform of <f)(x) at the magnetospheric boundary. The spatial variables 
x and z are dimensionless and run between 0 and 2?t. The drrift parameter is 8 = (v d L z )/(v A L x ), 
where L z and L* are the physical dimensions of the box. The magnetic flux function is, 


\j/(x,z) = T £ 




k =— oo ^1 + X. 2 & 


2Y 


172 exp 


j/Jfc jc±8( 


1 + X 2 * 2 \ 


1/2 


(26) 


The presence of the ionosphere means that the Alfven wave will be reflected. The 
solution for the case of a perfectly conducting ionosphere (<(> = 0) and a nonreflecting generator 
at the magnetospheric boundary is given by 


y(x,z) 


^ 1 ~ exp + 2 tc 5^1 + A. 2 A: 2 ) / Jjcos|^:z8(l + A. 2 £ 2 ) / 


( 27 ) 


where is the Fourier transform in x of a time independent arbitrary boundary profile f(x). 
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The steady state parallel electric field, the x-component electric field and the parallel 
cuirent density can be expressed in terms of the solution, Eq. (27). These are, 

E M (x, z) = -5 X 2 dl \|/ (x, z) = — 9^0 + 53 x \|/ (28) 

which is found from Eq. (10), the x-component of E is obtained from Eq. (8) and (22), 

E x U,z) = -8 _1 a 2 \|i(x,z) (29) 

and the parallel current density is given by Eq. (7) 

/,i(x,z) = 3 2 \|r(x,z) (30) 

There are several important and interesting features of this solution. Firstly it is obvious 
that the parallel electric field vanishes when either 5 = 0 or X, = 0. This means that for a 
significant parallel electric field to exist it is required that the arc be drifting normal to the 
direction of its alignment and that the transverse scale be of the order of the skin depth. The 
physical origin of the parallel electric field may be explained as follows. The quasineutrality 
condition requires 3 X J X + 9 Z J Z = 0, so that the cross field ion polarization current is balanced by 
the parallel electron current so as to maintain approximate charge balance. If, however, the 
electrons have finite inertia the parallel current they carry will lag slightly behind the cross field 
current, producing a small charge imbalance and thus a parallel electric field. Clearly the 
existence of the parallel electric field is implied by the fact that the region carrying parallel 
current is drifting into a current free region, hence for finite mass electrons an accelerating 
parallel electric field is required to maintain a steady state current in the frame moving with vj. 

Secondly, it is less obvious that E x vanishes if both 6 = 0 and X, = 0, but it does not 
vanish if either 5 or X is nonzero. This is physically reasonable, since only if there is no 
dispersion and no drift will the reflected electric field of the Alfven wave exactly cancel the 
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incident field so that the total field is zero. 

The steady state solution given by Eqs. (27) - (30) is intended to show the inter- 
relationship of the current and the electric field and it is probably a realistic description for dis- 
tances far from the ionosphere. An inhomogenous ionosphere, which was not considered in the 
analytic solution, significantly modifies the field structure near and within the ionosphere. The 
analytic solution given by Eq. (27) has been verified by comparing it to a two-dimensional time 
dependent simulation. The results are virtually identical. In the next section the three- 
dimensional time dependent model is extended to include a realistic inhomogeneous ionosphere, 
and numerical solutions are presented for a variety of situations. 

4. Numerical Simulations 

4.1 Inhomogenous Ionosphere Model 

The model equations (11) and (12) are expected to be valid for a homogeneous non- 
dissipative magnetosphere. A realistic ionosphere is both inhomogeneous and collisional; thus 
to extend the validity of the model into the ionosphere, additional terms must be incorporated to 
account for the effects of reflection from the density gradient, current closure via Pederson 
currents and resistive dissipation. To include transport terms, a coordinate transformation along 
the magnetic field is performed in order to resolve equilibrium ionospheric density variations. 

Assume n G (z) depends only upon z and take B 0 to be constant. Then define a new 
coordinate along B 0 through dz/dq = v A (z) for qj < q < q Q . Defining a(q) = [n 0 (q 0 )/n 0 (q)] 1 /2, 
the final equations are written in the frame moving with the generator are, 

[3, - v^a, + z x V<(> • V + v(^)]V 2 <), = [B c d q + a(q)z x Vy • V] V 2 y (3 l > 

(3, - v d d x + z x V(J> • V)(y - A,Vy) = V (32) 
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The ion-neutral collision frequency is v(q), and the resistive diffusion coefficient is T|(q). In 
this form the equations are suitable for an inhomogeneous ionosphere and they can be stably 
advanced in time. Without the re-scaling of the coordinate along the magnetic field, structures 
of a certain size along B will shorten relative to the spatial grid as they propagate into the region 
where the Alfven phase velocity decreases. This causes numerical instability due to fine scale 
structure accumulating at the grid scale size. 

For the purpose of acheiving Alfven wave reflection and current closure, very simple 
profiles of n Q (q), v(q), and n(q) are chosen. All time dependent simulation runs use 48 spatial 
grid points in the z direction. Let j denote the j th grid point in z starting at the bottom of the 
ionosphere. Then the profile is given by n 0 j = (256/ j2)3/2 , Vj = 50(16 - j)/15 , rjj = 0.2(16 - 
j)/15 for j=l,16; and n 0 j = 1 , vj = 0 , r|j = 0, for j=16,48. While not a completely realistic 
representation of the ionosphere, this simple profile accomplishes two things. First, Alfven 
waves are efficiently reflected due to the density gradient, and secondly the profile of the ion- 
neutral collision frequency is adequate to close most of the parallel current. The remaining 
parallel current is closed by applying the <J) = 0 boundary condition at the first grid point. It is 
possible to do a better job of modeling the ionosphere, however, for the purposes of the present 
paper and to do so would require many more grid points in the z direction which would make the 
three-dimensional simulations too expensive. In any case it has been found through 
experimentation with various other ionospheric profiles that the essential results are not changed 
as long as the criteria of good Alfven wave reflection and good current closure are met. 

42 The Numerical Code 

The three dimensional equations are solved using a mixed spectral and finite difference 
method with a second order Runge-Kutta as the time advance. The nonlinear terms in the 
transverse plane are computed using fast Fourier transforms and de-aliasing is implemented by 
grid shifting. The use of second order Runge-Kutta allows one to de-alias without additional 
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computational expense. The size of the x-y grid in all runs shown is 32 x 32, although some 
runs were performed using a 64 x 64 grid to determine if there were any qualitative differences. 
There were no significant differences. The z coordinate is finite differenced with the dependent 
variables <J) and \| f defined on alternate grid points to avoid instability associated with the decou- 
pling of the grid points. For conserving boundary conditions , the algorithm preserves the two 
global invariants discussed in section 2 to to better than 0.01% for 1000 time steps if the parallel 
and transverse Alfven Courant conditions are satisfied by a factor of two or more. 

The equations are solved in an appropriate set of dimensionless variables in which time 
is normalized to the parallel Alfven transit time. The total length of the system along the geo- 
magnetic field in code units is 2k. So that one complete transit along B takes 2k code time 
units. There is no fixed value of physical length along B, but the ratio of the x and z scale 
lengths is the same as the ratio of the transverse to parallel Alfven speeds and physically it is of 
the order of 1/1000. For the purpose of interpreting the simulation data, it is reasonable to take 
the the x -dimension to be about 10 km, the y dimension 20 km and the z dimension 10,000 km. 
To convert code time units to seconds, take the average Alfven speed to be 10 m/s, then 27t 
code time units is one second real time. These numbers should only be considered as guidelines 
since one is free to choose any set of scales, since the only physical parameter that selects a spa- 
tial scale is the ambient magnetosphere plasma density. For example, choosing the x-dimension 
to be 10 km and setting X = 1 in code units, implies the plasma density is about 10 cm' . Alter- 
natively if we take the plasma density to be 1 cm' J with X = 1, the the x-dimension of the simu- 
lation region is about 33 km. 

43 Steady state solutions 

All runs are initialized with zero field within the simulation box. The Alfven 
wave is launched into the box by application of the non-reflecting magnetospheric boundary 
condition Eq. (19), which is turned on by ramping the function f linearly in time until t = 2 tc at 
which time it is held at a constant amplitude of 2 for the duration of the run. 
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Contour plots of the electrostatic potential, the magnetic flux function, parallel electric 
field, and the parallel current density are presented in Figure 1 and 2 (clockwise starting at the 
top left panel) for the two-dimensional version of Eqs. (31) and (32) with 8 = 0.1, ^.= 1.0 for 
the S and V shock cases respectively. All plots are shown in the source drift frame. The results 
are for time t = 20, where t = 27t is the Alfven transit time along B. For times beyond t = 16 or 
so the results do not change significantly and thus a steady state has been achieved. 

There are several aspects of two-dimensional arc structure that can be seen from the con- 
tour plots. Regarding the closure of potential and parallel current contours, note that most of 
the potential contours close near the top of the ionospheric region which is consistent with 
electric field reflections from the density gradient. There are of, course, some potential lines 
which extend into the closure region (not drawn) to produce the electric field which drives the 
closure current. Most of the current contours, however, close well within the ionosphere, since 
the closure Pederson currents require ion -neutral collisions and higher plasma density. 

A comparison of the structure of the parallel current density and the parallel electric elec- 
tric field shows that the positive parallel electric field occurs in the region of negative 9 X J||. 
Thus the rightward drifting positive parallel electric field accelerates electrons downward to pro- 
duce the upward current which lags behind the parallel electric field. 

4.4 Three-dimensional time dependent simulations 

The two-dimensional steady state solutions do not remain two-dimensional if they are 
perturbed in the y coordinate. The results of three-dimensional nonlinear simulations are pre- 
sented for the same V shock and S shock magnetospheric boundary conditions as the two-di- 
mensional runs except that a small symmetry breaking perturbation is applied in the form given 
by Eq. (20). In order to see the effect of the source drift on the three-dimensional evolution, the 
drift parameter 8 is taken to be 0.0 and 0.1 for each of the V and S shock cases. A total of four 
runs are presented with the simulation parameters parameters given in table 1. 

In Figures 3 - 6 are the main results from the three-dimensional simulations correspond- 
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ing to runs 1 through 4 respectively. The differences between the non-drifting (S = 0) and 
drifting (5 = 0.1) runs are discussed first. The main qualitative difference is the asymmetry 
about the x = 7 t plane which is also present in the two-dimensional simulations. However, in the 
three-dimensional simulations for the S shock case, the asymmetry is also reflected in the 
qualitative differences between the negative and positive values of the potential and the current 
density in the x-y plane. For the V shock case, the main difference between the 6 = 0 run and 
the 5 = 0.1 run is that the non-drifting case has somewhat more structure than the drifting case. 
It was observed in other runs as well (not shown) that increasing 8 decreased the growth rate of 
the instability for both S and V shock boundary conditions. 

The structure along the magnetic field is much more field aligned than one would be led 
to believe from the x-z contour plots in the lower panels of Figures 3-6. This is because the 
total magnetic field, geomagnetic plus that due to Jy, has field lines that make a large angle with 
respect to the vertical in the y-z plane due to the strongly contracted z-dimension. Thus struc- 
ture existing in the x-y plane is projected onto the x-z plane. This is clearly seen by comparing 
the current density plot in the x-y and the x-z planes. Note also that there is much structure in 
the parallel electric field. The reason for this is the parallel electric field has one more derivative 
and hence more structure than the current density (see Eq. (28)). 

The source of the instability which is the origin of the three-dimensional structuring in all 
runs is the lowering of the magnetic energy through magnetic reconnection, i.e. a tearing mode. 
A tearing mode is implied for two reasons. One is because the local magnetic energy density 
exceeds the local flow kinetic energy and the form of the magnetic flux function and the flow 
stream function are similar. If the reverse were true then a shear flow instability like Kelvin- 
Helmholtz would be more likely. A stronger indication is that the hallmark of a tearing mode is 
found in the form of the current density viewed in the x-y plane. This is the formation of current 
filaments at magnetic x-points and flow vortices acting such as to force fluid into the magnetic 
islands between the filaments. The collisionless tearing instability as a possible source of 
auroral structuring was first suggested by Seyler , 1988. 
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In Figure 7 is a time sequence of contour plots of the electrostatic potential. The upper 
left plot is at time 16, the upper right at time 17, the lower left at time 18 and the lower right is at 
time 19. The entire sequence occurs on a real time scale of the order of a second. This is of the 
same order as that observed by Hallinan and Davis, 1970, and it is consistent with the estimated 
value of the tearing mode growth rate given in Sec. 2.3. 

The spectral characteristics of the instability in the nonlinear turbulent regime are quite 

interesting. The one dimensional power spectra of the potential (stream function) are presented 

in Figures 8-11. The first panel in each Figure is the x-direction power spectrum which is a 

function of k and an average over k . The second panel are the corresponding y-direction 
x y 

power spectra. Figures 8 and 9, which correspond to the S shock cases, are instantaneous power 
spectra at time 20, while Figures 10 and 11 corresponding to the V shock case are taken at time 
18. A qualitative description of the spectral behavior of the electric field spectral index, which is 
2 lower than the potential spectral index, may be stated as follows. During early stages of the 
instability, the spectra are much steeper (spectral index -3 to -5) than at later times when the 
flow can be considered to be more turbulent. At late times the spectral index very closely ap- 
proaches -5/3 , i.e. Kolmogorov scaling. The spectral data presented in Figure 8-11, corre- 
spond to the same times as in Figures 3 - 6. A comparison of the results of the S shock cases to 
the V shock cases shows that the spectrum is considerably steeper for the S shock cases. One 
must be careful in interpreting this observation, because spectra taken at later times for the S 
shock case (not shown) also shows a tendency towards a -5/3 power law. The correct 
interpretation of this result is that all cases approach a universal power law of -5/3, and the S 
shock cases shown represent spectra taken in an earlier stage of evolution. The S shock case is 
more slowly evolving because the amplitude of the magnetospheric boundary condition is such 
that the shear in the flow and the strength of the parallel current density is less than the V shock 
case, compare the values of the current density labels in Figures 1 and 2. 

The magnetic field spectrum has also been examined. It appears from the numerical re- 
sults that the linear relation Eq. (26) holds between the electrostatic potential and the magnetic 
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flux function; so that the electric amd magnetic field spectra are related as 



( 33 ) 


5 . Summary and Discussion 


5.1 Summary 

The primary conclusions that can be drawn from the analytical and numerical studies of 
the model are as follows: 

1. Oblique Alfven waves can form two dimensional steady state auroral arcs. For a 
nonzero collisionless skin depth there is associated a parallel electric field that drives the parallel 
current of the oblique Alfven wave. 

2. Density gradients as found in the topside ionosphere act as a good reflector of Alfven 
waves with an electric field reflection coefficient of T = -1. If the generator source is drifting 
relative to the plasma (oblique Alfven wave) then ionospheric reflections do not completely 
cancel and lead to a net transverse electric field. 

3. Two dimensional steady state arcs are unstable to three-dimensional perturbations. 
The arcs evolve through a combination of collisionless tearing and shear flow, with stronger 
driving field aligned current resulting in more rapid three-dimensional evolution. 

4. The transverse electric field spectrum of a three-dimensionally evolving arc appears 
to approach a k‘5/3 power law as time progresses. The approach is faster for higher shear and 
parallel current. At earlier times the spectrum may be quite steep. 

The results of the equilibrium model taken together with the results of the simulations 
indicate that the ionosphere play at least two principal roles in magnetosphere-ionosphere 
coupling. The first, and perhaps under appreciated, is the density gradients at high altitudes that 
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demark the upper boundary of the ionosphere (350-1000 km) can reflect Alfven waves quite 
efficiently. The combined effect of the incident and reflected oblique Alfven waves together 
with their response to electron inertia produces the electric field and parallel current structure of 
auroral arcs. The second (traditional) role is E-region current closure which closes the current 
carried by the electron beam that constitutes most of the upward parallel current. The electron 
current penetrates to much lower altitude than where Alfven waves are reflected and therefore 
has approximately zero parallel gradient, and hence no associated cross field closure currents. 
At such altitudes ( 110-160 km) where the ion-neutral collision frequency is sufficiently large to 
allow ion cross field currents is the circuit completed. 

5.2 Discussion 

The stationary arc solutions discussed in Secs. 3 and 4.3 clearly demonstrate the exist- 
ence of the parallel electric field, an estimate of the magnitude of E|| may be obtained from Eqs. 
(10) or (28). 


„ m e v d 

£|l - -^tMi (34) 

n 0 e 

3 16 3 

The quantities appearing in Eq. (34) are estimated as follows. 9 X ~ 2 x 10 m , n 0 ~ 10 m , 
J N ~ 10‘ 6 A/m 2 and guess the drift velocity to be v d ~ 3, 000 m/s in the magnetosphere which 
maps to a much lower velocity in the ionosphere. With these values, E|| ~ 0.2 mV/m, which is 
equivalent to one keV per 5000 km along B. The esimates of the parameters are typical, but any 
one may vary by a factor of two or more. It would seem to be possible to achieve a parallel elec- 
tric field of 1 mV/m and thereby accelerate electrons to a maximum energy of about 20 keV 
without requiring unrealistic magnetospheric conditions. 

Some comments are in order regarding observed electron precipitation data and possible 
connections with the model discussed here. The sounding rocket data of Arnoldy, 1970, 
showed rapid temporal variations of electron current and energy. This may have been be due to 
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internal reflections of the Alfven waves due to the magnetospheric density structure, Lysak, 
1988. It is possible to simulate auroral flicker using the model in this paper, but it would require 
a considerably more detailed magnetospheric density and magnetic field profile than was used. 

Recent data high resolution sounding rocket data, given in Earle [1988] and provided by 
C. Kletzing and R. Torben [unpublished], clearly shows the breakup arcs as inverted V’s of 
about 1 keV energy and a width of less than one kilometer. The electron precipitation satellite 
data of Gurnett and Frank, 1973, is not of sufficient resolution to detect the structure of the 
breakup arcs, as is the case for all satellite data that the author is aware. 

Whalen and Daly, 1979, have argued that static parallel electric fields will not produce 
the flat parallel electron distributions that they observed. In the present model the parallel elec- 
tric field is not static relative to the background plasma, but rather it drifts into a region of 
stationary electrons. So that electrons originating from different altitudes along the same mag- 
netic field line will acquire different energies by the time they reach the observation point below. 
Therefore, it does not seem that the data of Whalen and Daly are inconsistent with the model 
proposed here for electron acceleration. This point cannot be argued much further without a 
two-dimensional test particle simulation of electrons in the fields produced by a stationary arc. 
This topic is currently being investigated. 

The structuring found in the three-dimensional simulations bears resemblance to the ob- 
servations of Hallinan and Davis, 1970; however is difficult to assess just what simulation data 
should be compared to the observations of auroral curls. The potential contours bear the closest 
resemblance to the curls, but one might expect that the current density would be more appropri- 
ate for comparison. This is a rather difficult issue to resolve, since this kind of qualitative com- 
parison is highly subjective. Nevertheless it is fair to say that the collisionless tearing mode 
does have a structural morphology similar to auroral curls, as does the Kelvin-Helmholtz insta- 
bility, Keskinen et al, 1988, although the author does not believe that a shear flow instability is 
the origin of auroral curls in small scale arcs. The Kelvin-Helmholtz instability probably does 
play a role in large scale structuring associated with the low latitude boundary layer [Lotko and 
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Schultz, 1988]and the magnetopause[M/ura, 1984; Wu, 1986]. 

The spectral features of the electric field produced by the tearing mode process in the 
nonlinear regime appear to be universal taking on the form of a -5/3 power law. The DE 2 data 
oiBasu et al, 1988 upon a first consideration seems to be in only partial agreement with the sim- 
ulations. Three is, however, an interpretation that reconciles the apparent differences between 
the simulations and the DE 2 data as interpretated Basu et al. First a summary of the observa- 
tional results as characterized by Basu et al. They found for arcs having strong shear, the spec- 
tral index of the electric field was approximately -1.7, which agrees with the simulations; while 
for arcs having what was characterized as moderate shear, the spectrum was steeper in the range 
-2.5 to -3.5. A reconciliation of the differences with the moderate shear cases, is as follows. 
From the simulations it was found that lower shear and parallel current resulted in slower evolu- 
tion of the instability and steeper spectra at the same absolute time (recall the S shock spectral 
data). Therefore it is likely that the moderate shear spectra of Basu et al represents the instabili- 
ty in an early stage of evolution when the spectrum is steeper, while the strong shear arcs rapidly 
evolve towards the universal -5/3 spectral slope. Basu et al, 1988 also present density power 
spectra, which we cannot comment upon, since the present model does not have density evolu- 
tion. This is a subject for future research. 

The macroscopic fluid model described here has obvious limitations. Perhaps the most 
important of these is the simplified description of the electrons in which kinetic effects are mod- 
eled through the parameter X. Ions as well are treated as a low frequency fluid, so the physics of 
ion cyclotron waves is completely suppressed. Weak double layers (WDL’s) require a full ki- 
netic description of both species are are not described by the present model. It has been 
conjectured that WDL’s may be responsible for auroral electron acceleration, and there is some 
circumstantial evidence that the cumulative parallel potential drop of many WDL s may suffice 
to accelerate electrons to keV energy [Hudson et al, 1983; Koskinen et al, 1988; Tetreault, 
1988]. It is, however, the author’s belief that WDL’s are the result of electron drift and not the 
origin of the drift. The instability which forms the WDL’s (ion-acoustic or possibly something 
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else) may act to redistribute the field aligned potential drop into discrete packets, whereas the ac- 
tual origin of the parallel electric field is due to the oblique inertial Alfven wave. 
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TABLE 1. Parameters for three-dimensional simulations 


Run # 

1 S or V shock 

1 

1 

5 

1 

1 

A | 

X 

1 

1 

At 

1 

1 

1 S 

1 

1 

1 

0.0 

1 

1 

1 

17.8 . 1 

1.0 

1 

1 

1 

0.005 

2 

1 

1 s 

1 

1 

1 

0.1 

1 

1 

17.8 1 

1.0 

1 

1 

1 

0.005 

3 

1 

1 V 

1 

1 

1 

0.0 

1 

1 

1 

17.8 1 

1.0 

1 

1 

1 

0.005 

4 

1 

1 V 

1 

1 

0.1 

1 

1 

17.8 1 

1.0 

1 

1 

0.005 
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Fig. 1 Contour plots in the x-z plane of the electrostatic po- 
tential, magnetic flux function, parallel electric field and 
parallel current density for the S shock steady state arc, 
starting at upper right and going clockwise. 
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Fig. 2 Contour plots in the x-z plane of the electrostatic po- 
tential, magnetic flux function, parallel electric field and 
parallel current density for the V shock steady state arc, 
starting at upper right and going clockwise. 








Fig. 3 Contour plots of the electrostatic potential in the z=n 
plane, parallel current density in the z=7t plane, parallel 
electric field in the y = 7t plane and parallel current density 
in the y = 7t plane for run 1 starting at upper right and going 
clockwise. 






Fig. 4 Contour plots of the electrostatic potential in the z=n 
plane, parallel current density in the z=Jt plane, parallel 
electric field in the y = 7t plane and parallel current density 
in the y = k plane for run 2 starting at upper right and going 
clockwise. 






Fig. 5 Contour plots of the electrostatic potential in the z=rc 
plane, parallel current density in the z=7t plane, parallel 
electric field in the y = 7t plane and parallel current density 
in the y = k plane for run 3 starting at upper right and going 
clockwise. 









Fig. 6 Contour plots of the electrostatic potential in the z=7t 
plane, parallel current density in the z=k plane, parallel 
electric field in the y = k plane and parallel current density 
in the y = 7t plane for run 4 starting at upper right and going 
clockwise. 







Fig. 7 Time sequence of contour plots of the electrostatic 
potential in the y = 7t plane for run 3. The times are: upper 
left t = 16, upper right t = 17, lower left t = 18 and lower 
right t =19. 
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Fig 8 One dimensional potential power spectra for run 1 at 
time t = 20. Upper panel is in the x-direction and lower 
panel is in the y-direction. 
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Fig 9 One dimensional potential power spectra for run 2 at 
time t = 20. Upper panel is in the x-direction and lower 
panel is in the y-direction. 




log |potential|**2 log |potentiaI|**2 


Potential Spectrum vs kx 



Potential Spectrum vs ky 



Fig 10 One dimensional potential power spectra for run 3 at 
time t = 16. Upper panel is in the x-direction and lower 
panel is in the y-direction. 
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Fig 1 1 One dimensional potential power spectra for run 4 at 
time t = 16. Upper panel is in the x-direction and lower 
panel is in the y-direction. 





